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We study the properties of the eigenvalues of real random matrices and their products. It is known 
that when the matrix elements are Gaussian-distributed independent random variables, the fraction 
of real eigenvalues tends to unity as the number of matrices in the product increases. Here we present 
numerical evidence that this phenomenon is robust with respect to the probability distribution of 
matrix elements, and is therefore a general property that merits detailed investigation. Since the 
elements of the product matrix are no longer distributed as those of the single matrix nor they remain 
independent random variables, we study the role of these two factors in detail. We study numerically 
the properties of the Hadamard (or Schur) product of matrices and also the product of matrices 
whose entries are independent but have the same marginal distribution as that of normal products 
of matrices, and find that under repeated multiplication, the probability of all eigenvalues to be 
real increases in both cases, but saturates to a constant below unity showing that the correlations 
amongst the matrix elements are responsible for the approach to one. To investigate the role of the 
non-normal nature of the probability distributions, we present a thorough analytical treatment of the 
2x2 single matrix for several standard distributions. Within the class of smooth distributions with 
zero mean and finite variance, our results indicate that the Gaussian distribution has the maximum 
probability of real eigenvalues, but the Cauchy distribution characterised by infinite variance is 
found to have a larger probability of real eigenvalues than the normal. We also find that for the 
two-dimensional single matrices, the probability of real eigenvalues lies in the range [5/8, 7/8]. 
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I. INTRODUCTION 


The problem of the number of real roots of algebraic equations has a long history, but continues to attract attention, 
from the early works of Littlewood and Offord [lj to more recent developments [ 2,0 (see the latter article for more 
related history and references). Mark Kac, in a seminal paper, proved that the expected number of real zeros of 
polynomials (of order N) whose coefficients are chosen from a normal distribution of zero mean is ~ In IV Q. Others 
substantially extended these results and showed universality in that this is the leading order behavior, irrespective 
of the underlying distribution as long as they are zero-centered and have a finite variance |H|. Logan and Shepp 
0] studied the same when the coefficients are Cauchy distributed (and hence have infinite variance, and undefined 
average) the leading behavior goes as clnlV, with c ss 0.7413 which is larger than 2/tt k, 0.6366, and hence there 
are more real zeros in the Cauchy case than when the variance is finite. Thus, the number of real roots of a random 
polynomial are generally quite small. 

More recently, several works have explored the fraction of real eigenvalues for n x n matrices drawn from the 
real Ginibre ensemble. This ensemble consists of matrices whose elements are independently drawn from a normal 
distribution such as N(0, 1). For such matrices, it has been shown analytically that the expected number of real 
eigenvalues E n and the probability that all eigenvalues are real P„ j71 are given by 0-i 


,. E n . _ 
Inn — = \ 

n->oo y/n V 7r 


p _ 9 — n(n— 1)/4 

T7..7J. ^ 


( 1 ) 


Thus the probability that all eigenvalues are real tends to zero as the matrix dimension increases, although the 
expected number of them increases algebraically. 

Products of random matrices have also been studied, at least, since the classic work of Furstenberg and Kesten 
0 (see, 0,0 for a discussion of subsequent work and some applications). The study of the spectra and singular 
values of the products of random matrices is currently a very active area of research [13|, |l4[ , and recent progress has 
been reviewed in 0 ] . When studying a problem related to the measure of “optimally entangled” states of two qubits 
0 , the problem of the number of real eigenvalues of a product of matrices came up. It was shown in 0 that the 
probability of all eigenvalues being real for the product of two, 2x2 matrices from the real Ginibre ensemble is 7t/4. 
Somewhat surprisingly then there is a lesser probability that all eigenvalues are real for a single matrix (l/\/2, from 


Eq. (fljl ) than for a product of two. This observation motivated numerical explorations in 17] which showed that 


this probability tends to 1 as the number of matrices in the product increases. Numerical results showing that this 
is true for higher dimensions was also discussed therein. This result has been proven analytically by Forrester [181 . 
who calculated P njTl explicitly in terms of Meijer G-functions for square matrices, and generalised by Ipsen 0 to 
rectangular matrices. 

Most of the previous works on the properties of the eigenvalues have studied Gaussian matrices with independent, 
identically distributed (i.i.d.) entries. Exceptionally, Edelman, Kostlan and Shub [7] presented some numerical results 
indicating universality for the expected value E n of real eigenvalues ofnxn real matrices. However, to the best of 
our knowledge, there is no study exploring such extensions to the products of matrices. Thus, two conditions are 
relaxed: first, non-Gaussian i.i.d. elements are taken as entries of the matrices, and second, we consider products 
of such independent matrices. In this case, the entries of the product matrix are naturally correlated in a complex 
manner and their distribution is, in general, anyway non-Gaussian. 

Motivations for studying products of random matrices are well known in physics, and range from the study of 
localization in random media where transfer matrices are multiplied to dynamical systems where products of local 
stability matrices determine the Lyapunov exponents [0 . Applications for counting number of real eigenvalues of 
random matrices, especially of products, is of more recent provenance. The real eigenvalues of a class of random 
matrices indicate topologically protected level crossings at the Fermi energy in the bound states of a Josephson 
junction 0 |. While the relevant object is a single matrix and not a product, it is conceivable that this provides a 
context for further investigation. In the context of entanglement of two qubits [l7], the real eigenvalues of products of 
two matrices were relevant in finding the measure of optimally entangled states [161 ] . Different distributions of matrix 
elements would result in different sampling of states on the Hilbert space. Finally, it is worth pointing out that the 
fact that real eigenvalues dominate the spectrum of products of matrices implies that generic orbits of dynamical 
systems are not of the complex unstable variety. Complex instability can occur in Hamiltonian systems with more 
than two degrees of freedom [2lj . when exponential instability is combined with rotational action in phase space. 
However, as a consequence of discussions in this paper, it seems unlikely that they would generically arise for long 
orbits, as they depend on eigenvalues of products of Jacobian matrices being complex. 

To begin with, we ask: how do non-normal distributions affect the probability of real eigenvalues of a single matrix? 
In fact, a major part of this paper is concerned simply with 2x2 matrices and various families of the distributions 
of matrix elements, and the probability that the eigenvalues are real is calculated exactly in many cases in Sec. [TT] 
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In particular, we find the lower and upper bounds on the probability that both eigenvalues are real. In the family of 
symmetrized Gamma distributions, the range of probabilities is shown to be [5/8, 7/8]. The upper bound is associated 
with the case when the matrix elements have a large weight at the origin and the lower bound in the opposite case when 
the maximum is away from the origin. In fact, these same bounds are obtained in a class of truncated distributions 
thus indicating that the probability of real eigenvalues depends crucially on whether the distributions are concentrated 
near the origin or away from it. The lower bound of 5/8 is shown to be a tight one by proving this in general. The 
upper bound of 7/8 still remains specific to these distributions, however we do believe that this is widely applicable 
as well. 

We also find that the other features of the distributions of the matrix elements that play a significant role in 
determining the probability of real eigenvalues are the degree of smoothness and the existence of finite moments. It 
is noteworthy that these (apart from zero mean) are also the crucial features for universality in the case of random 
polynomials Given that the distributions are smooth and have finite moments, we tentatively propose that the 
normal distribution is the one with the maximum probability of finding real eigenvalues. If true, this provides a 
further unique characterization of the normal. 

Besides the bounds, we also observe a hierarchy for the fraction of real eigenvalues for matrices constructed from 
these distributions. Specifically, for some commonly occurring distributions, we find that the eigenvalues are more 
real when the matrix elements are chosen from Cauchy distribution than Laplace, Gaussian or uniform distributions, 
which are arranged in the increasing order of the decay. However, we must mention that such hierarchy patterns are 
not simply determined by the tail behavior of the distribution, and the trends appear to be more complex as explained 
in the following section. 

In Sec. M we numerically study the properties of the eigenvalues of the matrix obtained after taking product 
of several matrices, and find that the increase in the probability to one holds for several other distributions. In 
particular, we explore the zero-centered uniform distribution, the symmetric exponential (Laplace) distribution and 
Cauchy distribution, as simple representative ones. The hierarchy mentioned above is seen to hold even after the 
multiplication of independent random matrices. The expected number of real eigenvalues is also studied which tends 
to the matrix dimension exponentially fast (with the number of matrices in the product), although the rate of approach 
is distribution-dependent. 

It is unclear as to why the eigenvalues tend to become real when multiplying random structureless matrices. In 
an attempt at seeing how much this has to do with the act of multiplication, we study numerically the probability 
of real eigenvalues in the case of Hadamard (or Schur) products where the elements are simply the products of the 
corresponding elements of the multiplying matrices. Of course, in this case, the matrix elements remain uncorrelated. 
We find numerically that the probability of real eigenvalues increases with the number of matrices in the product, but 
tends to saturate at a value smaller than 1. Numerical evidence that this approach is a power law is also provided. 
This then highlights that the correlations built up in the process of (usual) matrix multiplication are responsible for 
the phenomenon that the fraction of real eigenvalues tends to one. 


II. REAL EIGENVALUES OF A SINGLE 2x2 MATRIX 

Let denote the probability that a product of K n x n random matrices has k real eigenvalues. In this section, 

we study the simplest case, viz., the probability that all the eigenvalues of a 2 x 2 matrix are real. We assume 
that the matrix elements are i.i.d. random variables chosen from the distribution p{x) with support on the interval 
[— u,u], where u is finite for bounded distributions and infinity for unbounded ones, and that the distribution p(x ) 
is symmetric about the origin as a result of which the mean of the probability distribution is guaranteed to be zero. 
We first give analytical results for the probability P^, which we study for many probability distributions, and then 

present some numerical results for the more general quantity P^J i n the following section. 

Consider a 2 x 2 matrix with i.i.d. elements defined as 

v x 

y w _ 

As the discriminant tr 2 — 4det must be nonnegative for real eigenvalues, the probability that all eigenvalues are real 
is 


pW _ 
r 2,2 ~ 


pU pU pU pu 


' —u J —u J —u j —u 


e[( v — w ) 2 + 4a :y\ p(v)p(w)p(x)p(y) dvdwdxdy. (2) 

Let the probability distribution oi z = v — w be q(z) with —2 u < z < 2 u. Observe then that the probability of real 
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eigenvalues can be written as 


P« _ L 

r 2,2 — 9 


pu p2u 

dx dy dz q(z)p(x)p(y)Q [z 2 + 4 xy] . 
1 Jo J 0 


( 3 ) 


Here the even symmetry of the distributions is used, and also that when the signs of both x and y are the same, 
the constraint imposed by the Heaviside function is trivially satisfied. In the above equation, the distribution q(z) is 
given by the convolution 


/ U 

p{v) p(w) S(z — (v — w)) dvdw 

-U 



p(v)p(v — z) dv. 


( 4 ) 


We note that this is indeed the usual convolution with z —> — z due to the symmetry of p(x). If z € [—2u,2u] the 
above expression for q(z) is valid, else it is zero. 

Equation (J3]) shows that the probability of both eigenvalues being real is at least one half for any probability 
distribution. The integral in Eq. © can be simplified if we use the convolution as 

^ pu pu p2u pu pu p 2 yjxy 

^ 2,2 = 9 + 4 / dx dy dz q(z)p(x)p(y) = 1 — 4/ dx dy dz q(z)p(x)p(y) . (5) 

z J 0 Jo J 2y/xy Jo Jo Jo 

The last equality follows from the fact that the convolution is itself a symmetric normalized distribution as 
q(z)dz = 1. Thus the evaluation of the probability of real eigenvalues reduces to evaluating the triple inte¬ 
gral above. 

It is also useful to write the distribution q(z) as 

-i POO pU 

q(z) = 2 nj dk ^ |iKA ° |2 ’ m = J P{X) ^ dX ' (6) 


Here p(k) is the characteristic function (or Fourier transform) of the probability distribution p(x). Next consider the 
integral defined as 


/ 0 pu p2u 

dx dy dz q(z)p(x)p(y)Q \z 2 + axy] . 
-u Jo Jo 


( 7 ) 


The value of interest is 1(4). Towards this end, differentiating 1(a) with respect to (w.r.t.) a, and performing the 
integral over the resulting delta function leads to 


dl_ 

da 


did 


pU Pli 

/ dx 
Jo Jo 


dy 


xy 


w 


2 y/axy 2 tt 


^p(x)p(y), 


( 8 ) 


where p(co) is given by Eq. ©. Noting that the above equation is valid only for a € (0,4), we integrate the last 
expression w.r.t. a from 0 to 4, and use 1(0 ) = 1/8 to get 1(4) and finally the probability that both the eigenvalues 
are real as 


P, 


(i) 


POO PU PU 

= 1—— / du dx dy\p(uj)\ 2 
7r Jo Jo Jo 


sin(2 ujy/xy) 
u> 


p(x)p(y). 


(9) 


Below we will apply the result in either Eq. ((5j or Eq. ([9|) to various choices of distribution p(x). Before turning to 
explicit calculations, we note that the unbounded distributions have the following scale invariance property. Consider 
a division of the random matrix elements by a nonzero constant b. If x' = x/b, then f) n , dx'p'(x') = 1 where 

p'(x') = bp(x) and v! = u/b. From Eq. Q, it is seen that the probability l is not affected by this scaling for 
unbounded functions, but for bounded ones, the limits should be redefined. While the above is an elementary analysis, 
we are not aware of it being discussed before. 

Under conditions of smoothness (in the sense that all the derivatives exist at least in intervals) and symmetry of 
p(x), it is interesting to enquire about the maximum value or upper bounds of P^-l , for instance. For the lower 
bound, we have already stated that the probability of all real eigenvalues is at least one half. A tighter lower bound 
is however 5/8, a proof of which was suggested to us by an anonymous referee which we discuss now. The probability 
of eigenvalues being real is 


P[(w — w) 2 + 4 xy > 0] 


-P[(u — w) 2 + Axy > 0 |xy > 0] 


w ) 2 + Axy > 0 |xy < 0], 


( 10 ) 
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Distribution po,„(x) ~ \x\ v 

Probability 

v = -4095/4096 

0.874959 

v = -7/8 

0.849868 

v = -1/2 

0.759836 

v = 0 

0.680556 

v = 1 

0.63709 

v = 3/2 

0.632888 

v = 2 

0.631023 

v = 3 

0.62928 

v = 4 

0.628361 

v = 200 

0.625078 

v = 400 

0.625039 


TABLE I. Probability that both eigenvalues are real for a 2 x 2 matrix with matrix elements chosen independently from the 
distribution in Eq. nm with [i = 0. As v approaches —1, the probability seems to limit to 7/8, while for v —> oo, it tends to 
5/8 (see text for details). 


where the factors of 1/2 arise as probability that xy > 0 or otherwise. However the first conditional probability is 1, 
and therefore 

P [(w — w) 2 + 4 xy > 0] = - + -P[(f — w) 2 — 4 xy > 0| xy > 0] (11) 

= - + ip[(u + w) 2 — 4vw — 4xy > 0| xy > 0, vw > 0] 

+ ^P[(u + w ) 2 — 4 vw — 4 xy > 0|a;y > 0, vw < 0] 

= - + yP[(u + w) 2 — 4vw — 4 xy > 0| xy > 0, vw > 0] 


( 12 ) 


2 4 

1 


H—P[(u + wf — 4vw — 4xy > 0|a;j/ > 0, vw < 0, \vw\ < xy] + 


1 


(13) 


4 

Here, in the second equality, further conditions on the sign of vw are used with equal probabilities for either occurrence 
(from the symmetry and independence of the distributions). In the third equality, a further certainty is carved out 
by using the possibility that when |vmj| > xy > 0 (occurence probability being 1/2), the discriminant is certainly 
positive. However our attempts at a similar approach for the upper bound were not successful. 

Although one can obtain get some insight into the bounds, for what kind of distributions these bounds actually 
occur is discussed below. 


A. Bounded distributions 


1. Symmetric Beta distribution 


Consider the symmetric Beta distribution with zero mean defined as 


PuA x ) 


r(/r + v + 2) 
2r(i + /z)r(i + i/) 


(1 - |x|r \x\ v , 


(15) 


with support on the interval [—1,1], This is a two-parameter symmetric family which can be nonsmooth at x = 0. For 
y = v = 0, the above distribution reduces to a uniform distribution, while y = 1, v = 0 and y = 0, v = 1 correspond 
to tent-shaped or V-shaped distribution respectively. 

Case: y = 0 

For the uniform distribution, the convolution q(z) = (2 — |z|)/4 on [—2,2] and zero elsewhere. Using Eq. (|5]>. we 
find that 


4$ O' = o) 


1 

2 



o 


V^ x v 


dz q{z ) = — = 0.680556. 


(16) 
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Distribution P/i,o(x) ~ (1 — |a:|) M 

Probability 

A* --1/2 

0.654534 

= 1/2 

0.695759 

CO 

II 

0.70085 

M = 1 

0.704854 


TABLE II. Probability that both eigenvalues are real for a 2 x 2 matrix with matrix elements chosen independently from the 
distribution in Eq. m3 with v = 0. The probability lies in the interval [5/8,11/15] for p € (—1, oo). 


Thus the uniform distribution results in a smaller fraction of real eigenvalues when compared to the normal distribution 
[jl, but not by very much. 

When p = 0 and v is nonzero and positive, the distribution po^(x) ~ |o;| I '0(l — |x|) is zero at the origin. As v tends 
to infinity, the weight of the distribution gets increasingly concentrated at ±1 and as Table [I] shows, the probability 
of real eigenvalues decreases. The distribution is smooth when v = 2k is an even integer. In this case, an analytical 
expression is possible for arbitrary fc, and given by Eq. (IA2I) of Appendix [A] Numerical analysis of Eq. (IA2I) for large 
v strongly suggests that 

P^^) = \+0(v~ 1 ) 1 i/»l. (17) 

An insight into the above result can be gained by the following heuristic arguments: If we consider the limiting 
distribution as the Bernoulli ensemble with only two choices of entries ±1 with equal probabilities, we get only 16 
matrices in the ensemble of which 12 have real eigenvalues. Thus it would seem that the ratio should have converged 
to 12/16 = 3/4. However, there are 4 cases where both the eigenvalues are exactly zero. If we believe that in the case 
of continuous distributions these are modified into cases with real and complex eigenvalues and equally, we get 6 cases 
of complex values and the fraction is consistent with 5/8 for real eigenvalues. Another related heuristic argument, 
closer to the continuous distributions under consideration, starts with the distribution 

p(x) = ^{S(x + 1) + 6(x-1)), (18) 

and applies Eq. © with p(co) = cos(w), leading to 

p , 1| = l_2 S 

* Jo u 8 

A proof of the assertion in Eq. o by analysing the double sums in Eq. (IA2I) seems rather difficult to obtain, but it is 
possible to tackle the Gamma distribution given by Eq. (El) , for which also the probability p(x) is zero at the origin 
and increases algebraically away from it, and show that indeed the limit probability is 5/8, see Sec. Ill B 21 Thus the 
lower bound derived above is a tight one. 

We now turn to the case when — 1 < v < 0 where the probability distribution piles up at the origin. As shown 
in Table Q] the probability of real eigenvalues increases as v decreases towards — 1. It is obvious that if the matrix 
elements are chosen from a Dirac-delta distribution centred about zero, both the eigenvalues are definitely real. But 
whether the probability limits to something less than one as v approaches —1 is of natural interest. As shown in 
Appendix [Al we find that 

P^M = l~0(l + v), u^-1. (19) 

Thus it is interesting that the distribution ~ \x\ u restricted to an interval spans a range of behaviors for the 
probability of real eigenvalues. As the matrix elements probability gets increasingly piled up at the ends (±1), the 
probability of real eigenvalues decreases and tends to 5/8, while in the opposite case when the elements are piled up 
around origin, the probability approaches the maximum value of 7/8. 

Case: v = 0 

For v = 0, the probability distribution p^. v { 0) ^ 0 and we find that with increasing p 7 the probability that the 
eigenvalues are real increases, refer Table [II] Partial analytical results are possible. For example, for the case when 
p = 1, v = 0 (“tent” distribution), the convolution is 


q(z) = 


i(2^M) 3 1<M<2 

1 ( 4 - 6|^| 2 + 3\z\ 3 ) 0 < 1^1 < 1 . 


( 20 ) 
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The resulting integral in Eq. ([Sj can now be done using hyperbolic coordinates. This results in 


P±% = l,u = 0) 


16143 23 In 2 

22400 ~~ 1008 


0.704854. 


( 21 ) 


It is reasonable to expect that as p increases (and v = 0), the probability increases to that when the elements are 
distributed according to the Laplace distribution ~ exp(— \x\) which is obtained as p —> oo. We will deal with the 
Laplace distribution below, but state here that the probability of real eigenvalues in this case is 11/15 ~ 0.733. 
However, for negative p, where the matrix elements have a tendency to be close to ±1, from the discussion in the last 
subsection, we expect the probability to approach 5/8 as p —> —1. Thus for the distribution p^^^x), the probability 
of real eigenvalues lies in the range [5/8,11/15]. 


2. A smooth family 


Consider another class of smooth zero mean distributions bounded on the interval [—1,1] defined as 


p v (x) 


rfo + §) 

+ !) 


(l-x 2 y , 77 > 1 . 


( 22 ) 


Of course, for 77 > 0, the distribution is continuous but not smooth at ±1, but this does not seem crucial. When 
77 = 0 , Prj{x) is uniform on the said interval, and as p — > 00 the distribution approaches the normal distribution with 
the variance scaling as I/ 77 . As discussed above, the probability of real eigenvalues is independent of the variance 
and therefore, we expect that the large Rvalue for this probability will coincide with the known result for the normal 
distribution, namely, l/y/2 ~ 0.707 • • • [2J (also, see Sec. Ill B II) . The probability of real eigenvalues is calculated for 
some values of integer 77 in Appendix [Bj and we find that the probability indeed approaches l/\/2 with increasing 77 . 

When p is negative, the distribution p{x) diverges at x = ±1, and from the discussion in the preceding subsection, 
we expect the probability of having real eigenvalues to approach 5/8 as 77 —f — 1. The case of 77 = —1/2 is that of the 
arcsine distribution. The convolution with itself can be found and is a complete Elliptic integral. However, even if 
further analytic results seem to be hard, this enables a more accurate numerical estimate of the probability of real 
eigenvalues which is ~ 0.662. Further decreasing p makes the numerical evaluation unstable, but results indicate a 
monotonic decrease in the probability. 


B. Distributions with infinite support 

1. Gaussian 

The case of normal or Gaussian distribution is the most studied and there are general results [7l.H7l.fl8l]. We consider 
a zero centered, unit variance Gaussian distribution for the matrix elements given by 

e -x 2 /2 

P( x ) = x € (- 00 , 00 ). (23) 

V Z7T 

The probability of both eigenvalues being real has already been shown, using Eq. ©, to be exactly equal to l /\/2 in 
E3. and alternative derivations naturally exist, see the works of Edelman 0 and Forrester [l8j- However to place it 
in the context of this paper, we rederive the probability of real eigenvalues in this case in Appendix [Cl 

It is well known that the normal or Gaussian distribution is singled out in numerous ways, for example, as one 
that maximizes entropy for a given mean and variance, or as the limit of sums of random variables. In the context of 
the present work, it seems plausible that the normal distribution is once again to be singled out as the distribution 
of matrix elements that maximizes the probability of finding real eigenvalues among the class of smooth, symmetric, 
and finite moments distributions. To test this further, we have looked at distributions p(x) such as ~ exp(—x 4 ), 
~ exp(— x 2 — rx 4 ) (r > 0 ) etc. and verified in all these cases that the probability of real eigenvalues is indeed less than 
l/\/2. We advance this proposition tentatively based on evidence gathered so far, and from results to be presented 
below. 








2. Gamma distribution 


We next consider the symmetrized Gamma distribution defined as 

Pi 0 ) = > 7 > °- 


(24) 


where the scale for the exponential decay has been set to one. The special case of Laplace distribution for which 7 = 1 
and the general case using the convolution route are discussed in Appendix [D] The convolution itself is given by 


q{z) = 


4T(2 7 ) 


e-HUI 2 ''- 1 - 


2^+i (7) 


\z\^K^_A\z\), 


(25) 


where K v (z) is the modified Bessel function of the second kind [ 22 }. The results obtained by analysing the resulting 
integrals are shown in Table EU and we find that with increasing 7 , as the weight of the distribution at the origin 
decreases, the probability P ^ also decreases. 

The behavior of the probability of real eigenvalues as 7 approaches zero can be understood as follows. As 7 
approaches zero from the positive side, the distribution and both the terms in the convolution q(z ) diverge at the 
origin as \/\z\. The latter can be seen as K_i/ 2 (z) = K-\ / 2 {z) = \fnj2e~ z / yfz, and lim^o T{2z)/T(z) = 1/2, as 
follows from the duplication formula for the Gamma function [2(31 . Since we expect that the divergence at the origin is 
all that matters, the probability for real eigenvalues will converge to the case of the bounded distribution \x\ u already 
studied above, and therefore the probability for real eigenvalues increases to 7/8 as 7 —► 0. 

For increasing 7 > 1, the distribution itself is peaked away from the origin and has two symmetric maxima at ~ 7 . 
The probability P 2 l 2 ( 7 ) now decreases from the 7 = 1 value and monotonically seems to approach 5/8, the lower limit 
for the bounded \x\ v distributions as v increased. Indeed the piling up of the probability of the matrix elements at 
two symmetric “walls” makes this plausible. To analyze this further, we restrict attention to integer values of 7 and 
give an exact evaluation in terms of finite sums as 


p 2 ( ,2 ( 7 ) = \ + ^ 1 ( 7 ) + ^ 2 ( 7 ), (26) 

where A\ and A 2 are given by Eq. (ID8I) and originate from the first and second terms of the convolution q(z ) in 
Eq. (1251) respectively. To find the limiting value for the probability of real eigenvalues for large 7 , we note that like 
the distribution p 7 (x ), the first term in the convolution q(z) is also peaked around z ~ ± 2 y, but the second term 
peaks at z = 0 so that the two terms in the convolution have practically disjoint supports. As a consequence, for large 
7 , the dominant contribution to the probability in Eq. (1261) comes from A\{p/) which represents the overlap between 
the distribution and the convolution, and ^ 2 ( 7 ) can be neglected. As discussed in Appendix [D] on analysing the sum 
A 1 ( 7 ), we get 


P StM = l + + 7»1. (27) 

In summary, for the symmetrized Gamma distributions, the probability of real eigenvalues also seems to be in the 
range [5/8, 7/8]. The distribution having a nonanalyticity at the origin leads to larger probability for real eigenvalues 
compared to the normal distribution and it seems comparable to the bounded power law distributions \x\ v studied 
above except in the rates of convergence. 


3. Power law distributions 

A qualitatively different process is interesting to consider, and as in the case of random polynomials 0, it will 
be interesting to study what happens when the underlying probability distributions have diverging moments. The 
Cauchy distribution which is the simplest and best studied of these and occurs in many contexts, is given by 

(28) 

for x G (— 00 , 00 ). Again the possible parameter in the distribution is rendered inoperative via scaling. Using 
p(w) = e ~and performing a change of variables x = u 2 in ([5jl. the integral w.r.t. u can be calculated. Then 
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Distribution symmetric Gamma 

Probability 

7 = 1/4 

0.824051 

7 = 1/2 

0.784155 

7 = 1 

0.733333 

7 = 2 

0.68325 

7 = 3 

0.660393 

7 = 10 

0.633238 

7 = 100 

0.627494 


TABLE III. Probability that both eigenvalues are real for a 2 x 2 matrix with matrix elements chosen independently from the 
symmetric Gamma distribution in Eq. (HU). The probability lies in the interval [5/8, 7/8] for 7 £ (0, 00 ), with the lower bound 
corresponding to 7 = 00 and upper to 7 = 0 . 


computing the integral w.r.t. ui using series expansion yields 


P. 


(i) 



1 

1 + y 2 


tan 1 



(29) 


Although the integral is written as if it is carried out, it is in fact a numerical evaluation which is almost certainly 
correct. Evaluations using the convolution path lead to interesting alternative integral forms of 3/4, but none of them 
(including the above) seem to be either in standard tables or calculable using symbolic mathematical packages. 

In general, one may consider the distribution p a {x) ~ 1/(1 + x 2a ) , a > 1 which possesses finite (and nonzero) 
moments up to order 2a — 2 only. Numerical evaluation of Eq. (0 gives 0.7076005 and 0.694185 for a = 2 and 3 
respectively, both of which are smaller than the result obtained above for the Cauchy distribution which is the slowest 
decaying power law distribution with finite mean. We also note that compared to the Gaussian case, the eigenvalues 
are less likely be real for the power law-distributed matrix elements with a > 2, and therefore the probability P^ 
is not merely determined by the decay behavior of the distribution of the matrix elements. Other distributions such 
as 1/(1 + |a;| a ) , 1 < a < 2 which is slower than Cauchy distribution but not smooth, or power laws corrections that 
vanish or diverge at the origin to the fat-tailed distributions have not been investigated. 

To summarize, the case of a single 2x2 matrix with elements drawn from various i.i.d. distributions have revealed 
an interesting phenomenology. It may seem like the weight of the distributions near the origin matters and if the 
values are clustered around zero, the probability of real eigenvalues increases. This statement must however be 
qualified: the normal distribution with however a small variance always have only 1 / y/2 probability of having real 
eigenvalues. Thus the differentiability of the underlying distributions play an important role. In the extreme case 
of delta distributions, we may have 100% eigenvalues real. But given that the distributions be smooth to all orders, 
the normal distribution seems to be singled out as the one with the largest probability of real eigenvalues. Also, the 
probability of real eigenvalues seems to be in the range [5/8, 7/8] for the classes of distributions considered here. For 
the case of the Cauchy distribution which is smooth but has diverging variance, the probability is large at 3/4 for 
real eigenvalues. Thus this rather simple problem of the probability of real eigenvalues of random real 2x2 matrices 
seems to possess a multitude of interesting features that warrants further study and clarification. 


III. REAL EIGENVALUES OF PRODUCT OF MATRICES 

We now turn to the properties of an n x n matrix obtained after taking a product of K square matrices, and study 
how the probability that some or all of the eigenvalues are real behaves for K > 1. It is of interest to see if the results 
for single 2x2 matrix described in the last section carry over to the product of matrices. Thus there is a two-fold 
generalization: (1) products of matrices are considered, and (2) their dimensionality can be more than two. 


A. Asymptotic value and maintenance of hierarchy 

We numerically studied the probability that k eigenvalues are real for a product of K n x n random matrices. 
As the products of random matrices can have, in general, a positive Lyapunov exponent 0, numerical procedures 
for all cases of products “renormalizes” the matrices for each product by dividing the Frobenius norm. This however 
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FIG. 1. Comparison of probability that all eigenvalues are real for a product of K random matrices with different symmetric 
distributions and the dimensionality n = 2 (main) and 8 (inset). The plot is based on 10 J independent realizations. 


does not alter the quantities of interest here. For most of the discussion, the matrix elements are assumed to be 
i.i.d. random variables distributed according to one of the following symmetric probability distributions: uniform, 
Gaussian, Laplace and Cauchy. Note that all of these distributions are finite at the origin, but have different tail 
behavior. The data were averaged over 10 5 — 10 6 independent realisations of the random matrices. For the 2x2 
matrix, there can be either zero or two real eigenvalues, but for higher dimensional matrix with n = 8 that we consider 
here, the probability of real eigenvalues is nonzero for k = 0, 2,4 ,6 and 8 only. Numerical results for the probability of 
all real eigenvalues for n = 2 and 8 are presented in Fig. 0 for various probability distributions as a function of the 
number of matrices in the product. We find that the probability of all eigenvalues being real increases to unity with 
K monotonically for most distributions (see, however, the case of Gamma distribution defined by ( 1 M 1 ) with 7 = 10 ). 
This effect has been previously observed by one of the authors G3 for Gaussian-distributed matrix elements. Here 
we find that this result is quite general in that the eigenvalues tend to become real with increasing K when matrix 
elements are distributed according to non-normal distributions as well. 

The second important point illustrated by Fig. 0 is that the same hierarchy as observed for the K = 1 case 
continues to hold for K > 1. Explicitly, the probability of all eigenvalues being real for the uniform distribution is the 
smallest followed by the Gaussian distribution, the Laplace distribution and finally the Cauchy distribution. Thus 
the slowly decaying distributions appear to have larger probability of having all real eigenvalues. This feature is seen 
to be valid for higher dimensions as well (results not presented). Thus the hierarchy apparent even with a single 2x2 
matrix continues to hold for larger dimensional matrices as well as the product of the random matrices. Furthermore, 
Fig. 0 shows that the ordering of the probability of real eigenvalues according to their tail behavior is not special 
to the probability of all real eigenvalues, but holds when k n as well. Note, however, while the probability of all 
real eigenvalues increases monotonically towards unity, the probability of k < n real eigenvalues decays to zero, as 
the number of matrices in the product increases (see, the inset of Fig. ©)• 

The increase in the probability of all real eigenvalues is reflected in the average number of real eigenvalues as well 
which is given by E^ = YJk-o • Since, as discussed above, the probability P^ K J —> 5 n ,k as K increases, 

the average E^ i ' > approaches the dimension of the matrix. Numerical results in Fig. © for 8 x 8 matrices strongly 
suggest an exponential approach of the expected number of real eigenvalues to the dimension of the matrix. The 
data also points to the continuance of the hierarchy even for the average number of real eigenvalues indicating a more 
microscopic adoption of the hierarchy to all . 
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FIG. 2. Probability that 4 eigenvalues are real for a product of K 8 x 8 random matrices whose elements are chosen from 
uniform, Gaussian, Laplace and Cauchy distributions, each with zero mean, based on 10 6 realizations. The inset shows the 
probability that k eigenvalues are real for a product of K 8 x 8 zero mean, Laplace-distributed random matrices. 


B. Effect of correlations between matrix elements 


1. Hadamard product 


We now consider the Hadamard (or Schur) product wherein the elements of the product are simply the products 
of the corresponding elements. Thus if o is used to denote the Hadamard product of two matrices A and .B, then 
(A o B)ij = AijBij. The Hadamard product of two positive semidefinite matrices is also positive semidefinite, unlike 
the usual product, and is one of the reasons it is important J24[. However the reason why we choose to study Hadamard 
products is that although the matrix elements of the product matrix are products of random numbers, there is no 
correlation amongst the matrix elements. Thus one can hope to disentangle two possible mechanisms that may be 
responsible for the phenomenon that all the eigenvalues of a random product matrix tend to be real: the act of simple 
multiplication which is presumably causing the matrix elements to have large weight near the origin and the various 
addition of such products that are leading to correlations between the matrix elements. 

To this end, we consider the eigenvalue properties of a Hadamard product of K random nxn matrices with elements 
distributed according to uniform, Laplace, Gaussian and Cauchy distributions, each with zero mean. Figure (J3J shows 
the case of Hadamard products of 2 x 2 matrices, and we observe that the hierarchy effect is maintained at any 
K (although there are more fluctuations in this case as compared to that of the ordinary matrix product). But, 
importantly, we find that the probability does not approach unity with increasing K. 

We can get some insight into this result by calculating the distribution of the matrix elements of the Hadamard 
product matrix for some distributions, and appealing to the results for the 2x2 matrices obtained earlier in Sec. HU 
Let us first consider the case of Hadamard product of matrices with elements drawn from the bounded distribution 
Po, v (x) defined in Eq. (fl5l) . Let zk = X\X2...xk represent a random variable formed by the product of K i.i.d. random 
variables. As discussed in Appendix [0 the distribution of the product for this case is given by 


Pk(zk = X\X2—Xk) 


{v + l) K \z K \ v 
2(K — 1)! 



K -1 

0(1-M). 


(30) 


For v < 0, the above distribution diverges at the origin while for positive v, it vanishes at \z\ = 0 and 1 and is a 
nonmonotonic function of \z\. Thus, except for the uniform case, the behavior of the distribution pk(zk) is similar 
to that of po,v near the origin. Using the above result, it is possible to numerically evaluate the integral in (J5]) to 
obtain for various values of v and K. For example, for v = 0, we find the probability of real eigenvalues to be 
0.738779, 0.767331, 0.782558, 0.792032, 0.798561, 0.803376 for K = 2,3,4, 5, 6, 7 respectively. For larger K, numerical 
integration does not converge; however, as shown in Fig. [U for some representative values of z/, the data obtained 
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FIG. 3. Expected number of real eigenvalues for a product of K random 8 x 8 matrices with elements chosen from uniform, 
Gaussian, Laplace and Cauchy distributions, each with zero mean. The plot is based on 10 6 realizations. 


from direct sampling indicates the approach to a probability less than 7/8. For the case of uniform distribution, 
the probability seems to saturate around 0.84. But for negative and positive v 1 the probability of real eigenvalues 
approaches a value higher and lower than 0.84 respectively. This pattern is consistent with the results in Sec. HD where 
the probability of real eigenvalues decreased with increasing v. 

We next consider unbounded distributions, and start with the case of Hadamard product of matrices with Gaussian 
distributed elements with zero mean and unit variance for which fl 8 l . |25| 


1 poo pOO 

P2{z 2 ) = — / / € r x ll 2 e r x ll 2 S(z 2 - Xix 2 ) dxidx 2 

J —oo J —oo 

4f r w-., il = « ae( .„,4 

2tt J_ 00 tt 


(31) 


which diverges logarithmically at the origin since Kq{z) —> ln(l /z) as z approaches zero. For this distribution, the 
integral in Eq. (0 can be numerically evaluated by using the fact that its convolution is given by a Laplace distribution 
(shown in Sec. Ill B 2) yielding P 22 = 0.757164, which is consistent with the value obtained in Fig. 0 by direct 
sampling. For K > 2, the probability distributions are given by the Meijer-G functions [l^, [25| making it difficult to 
handle this using Eq. ©. We therefore give the results obtained via numerical sampling in Fig. [4] which again indicates 
an approach towards 0.84. In the case of Hadamard product of two matrices with Laplace-distributed elements with 
zero mean given by (Phil) with 7 = 1 , we have [25| 


1 pOO poo 

P2{z 2 ) = - / e -l Xl l e -l X2 l S(z 2 - xix 2 ) dx\dx 2 

^ J —00 J —00 


1 

2 


dx 1 


K 0 (2y/\z^\), z 2 G (- 00 , 00 ), 


(32) 


where the last step in Eq. (1321) has been done by a series of two transformations, x\ —> x \followed by x\ —> e Xl . 
Numerical evaluation of Eq. 0 using this yields P 2 2 2 = 0.773849, and the results obtained from direct sampling 
show a convergence towards a value close to 0.84. For the Cauchy-distributed matrix elements, the distribution of 
the product of two random variables is known to be [26j 


P2(z 2 ) 


1 nzj 

^{4-lY 


(33) 


which diverges at the origin logarithmically. 
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Thus, as Fig. Q shows, the probability of all real eigenvalues monotonically increases with K and saturates to a 
value in the range 0.82 — 0.84. Note that these numbers are less than 7/8 which, from our analysis in the previous 
section, is expected of distributions that diverge as a power law, | x\ v with v —> —1, at the origin. From the above 
calculations for K = 2, it is clear that the uniform distribution and the three unbounded distributions discussed above 
have Hadamard product elements distributed according to a probability distribution that diverges as ln(l/|x|) at the 
origin. In fact, the elements of the Hadamard product of uniform distribution as well as Cauchy distributed random 
matrices have probability distributions that diverge as \n{l/\x\) K ~ l at arbitrary K [26j. It is reasonable to expect 
that the Gaussian and Laplace cases also diverge in a similar manner. For these cases, as the divergence at the origin 
is somewhat “slower” than a power law, it seems reasonable that the probability of eigenvalues being real is a shade 
smaller than that of single 2x2 matrices whose elements are distributed according to the power law discussed above. 

Assuming that the probability of Hadamard products saturate to the same value for these four distributions in the 
K —> oo limit, and taking this to be ss 0.846 (in the case of n = 2), our numerical results show a power law approach 
to the constant: 


p(K) 

r 2,2 


p(°°) 

r 2,2 


c_ 

K°' 


(34) 


where C is a positive constant and the exponent 8 m 0.675, 0.649, 0.654, 0.621 for uniform, Gaussian, Laplace and 
Cauchy distribution respectively (see the inset of Fig. ljl)l). Thus, here the probability of real eigenvalues approaches 
the asymptotic value as a power law in comparison to the ordinary matrix product case, where this approach is 
seen to be exponentially fast 0 - The exponential behavior is also seen in the expected values of real eigenvalues 
as in Fig. [3] We also looked at the fraction of expected number of real eigenvalues for higher dimensional matrices 
(n = 2,3,4) when each matrix in the Hadamard product has Gaussian-distributed elements with zero mean, and find 
that the fraction of real eigenvalues increases and saturates to a value less than one, unlike that observed for the case 
of usual matrix products. We also find that the asymptotic value of this fraction decreases with an increase in the 
dimensionality of the matrices. 

The convergence of the probability that all eigenvalues are real as K —> oo for the Hadamard product can also be 
understood using the Central Limit Theorem applied to the logarithm of the absolute value of the product random 
variable. The elements of the Hadamard product matrix would then admit a symmetrized log-normal probability 
distribution of the form 0] 


Pk(z K = XiX2--X K ) = — - 

21 z K 

for large K , with the parameters p and er being the mean and standard deviation of log|x|. Numerically, we find 
that the probability that all eigenvalues are real for a matrix with elements distributed according to Eq. (1351) indeed 
converges to a value around 0.84 at large K, as for the Hadamard product. 


\V%7 tKo^ 


exp 


(log|zftr| - Knf 


2 Ka 2 


, z K e (- 00 , 00 ), 


(35) 


2. Usual matrix product 

We now return to the properties of the matrix obtained after multiplying K matrices, and ask why the asymptotic 
probability of real eigenvalues saturates to unity unlike that for the Hadamard product. We can think of two possible 
reasons: one, the (usual) product matrix elements have probability distributions that are significantly different from 
the Hadamard case. The other possibility is that the correlations between the matrix elements in the case of usual 
products lead to an increase in probability of all eigenvalues real to unity as opposed to the Hadamard case where 
this is not true. It is not obvious whether the difference in the probability distribution or the correlations between 
the matrix elements is responsible for a higher P 22 ). 

To gain some insight into this intriguing question, we looked at the matrices with i.i.d. random variables as elements, 
each distributed according to the probability distribution of the product matrix elements at various K. Such matrices 
can be easily generated by taking independent samples of usual product matrices at each K and forming new matrices 
by picking the corresponding matrix elements from the independent product matrix samples. These new matrices 
would then have elements with probability distributions corresponding to that of elements of a usual product matrix 
at product length K } but the matrix elements would now be independent, having been selected from independent 
product matrix samples. Fig. (|5J) shows a comparison between P^P for the Hadamard product and usual product 

and P 22 f° r matrices with probability distribution of elements corresponding to that of usual product at product 
length K , but without any correlations between the elements. It is clear that the difference in probability distribution 
of matrix elements between the Hadamard and the usual product has an effect of lowering for the usual product 
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FIG. 4. Comparison of probability that all eigenvalues are real for Hadamard products of K 2 x 2 random matrices for some 
symmetric distributions based on 10 s realizations. For the Beta distribution with v — 3, the x-axis is scaled up by a factor 
5 since the convergence to the asymptotic value occurs for very large K. The inset shows the power law approach of the 
probability of all real eigenvalues to the asymptotic value which is less than unity, for the Gaussian case. The plot is based on 
10 5 realizations for the Beta distribution, and 10 6 for the rest. 


case to a value below that of the Hadamard product case. However, the presence of correlations between the matrix 
elements of the usual product matrix is seen to have an effect of a substantial increase in P 2 ^P, finally leading to an 
asymptotic value of one. 

The effect of correlations can be seen analytically for the case of Gaussian-distributed matrices when K = 2. We 
have seen that the distribution of the product of two Gaussian-distributed random variables with zero mean and 
unit variance is given by Eq. (13ID and the probability of real eigenvalues is 0.757164. The distribution of the matrix 
elements obtained on taking the usual matrix product corresponds to that of the sum of two random variables, each 
distributed according to the distribution of the product of two random variables. The characteristic function of 
Eq. (1311) is given by p(k) = (1 + fc 2 ) -1 / 2 , k £ [—1,1] (due to Eq. (11.4.14) of (111)- It is easily seen that this p(k) 
is the square root of the characteristic function of Laplace distribution. Hence, the product matrix elements have a 
Laplace distribution, for which we know from Eq. (ID2D that P 22 would have been 11/15 = 0.7333, had the elements 
been independent. But, the correlations between the matrix elements raise this probability to 7 t/ 4 = 0.7854 jl7l fl8j 
for the usual product, hence presenting a strong evidence that the correlations between the matrix elements are in 
fact leading to P 2 *P —I 1 observed for usual products. 

We also measured the correlations between the matrix elements to see how correlations increase with increasing K 
for 2x2 matrices with Gaussian-distributed matrix elements. More precisely, we consider 


- {xi)(xj) \^ 1/K 


i = 1, ...4 


(36) 


where Xi are the matrix elements obtained after taking the product of K matrices. For K = 1, obviously Ci = Si y i. 
But with increasing K , all the four Ci s seem to be of same order i.e. the variance of the new distribution and 
the correlations are similar. In particular, we find that for K = 2, the correlations Ci = 1.417,0.047,0.031,0.057 
for i = 1, 2, 3,4 respectively, but these numbers increase to 1.885,1.446,1.349,1.462 and 1.905,1.500,1.475,1.667 for 
K = 12 and 20 respectively. These data suggest that similar to the probability P 22 , the correlations also tend to 
saturate with increasing K. 
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FIG. 5. Comparison of P 2 2 ^ f° r Hadamard products and usual products with correlated elements and P 2 for matrices with 
i.i.d elements distributed according to probability distribution of product of K independent matrices, for 2x2 zero mean, 
Gaussian-distributed random matrices. The plot is based on 10 6 realizations for the Hadamard product, and 10 s for the rest. 


IV. DISCUSSION 


How does the probability of real eigenvalues for a 2 x 2 random matrix depend on the probability distribution of the 
matrix? It is quite surprising that such an apparently elementary question presents many novel challenges, and throws 
up some surprises. Here we find that the probability of real eigenvalues depends on several detailed features of the 
probability distribution, some of which we have identified here as the finiteness and smoothness at the origin, existence 
of maxima away from origin and the finitness of the moments. We have shown that eigenvalues are most likely to 
be real for distributions limiting to p(x) ~ l/|x| which is the most divergent distribution at the origin (with zero 
mean) and occur with a probability 7/8 = 0.875. That the large weight at the origin correspond to high probability 
is perhaps not surprising since when the matrix elements are distributed according to <5(a;), this probability is unity. 
However, interestingly, here we find the maximum probability to be less than one. The origin of this naturally rests 
in the smoothness of the probability distributions considered. We also mention that the values are not universal; for 
example, the probability turns out to be different for the Beta distribution po, 2 (x) and for Gamma distributed matrix 
elements with 7 = 3, although both have the same behavior near x = 0. 

Our results suggest that the probability that all eigenvalues are real is larger for distributions with large weight 
at the origin and that decay slowly. More concretely, there exists a hierarchy between the probability values of the 
different distributions, which is Cauchy (0.750) > Laplace (0.733) > Gaussian (0.707) > uniform (0.6805). Moreover, 
for class of distributions that decay in the same manner, the distributions with higher weight close to zero have a 
higher probability of real eigenvalues. For example, although the probability p^l for distribution Eq. (1241) for 7 < 2 
is larger than the uniform distribution as it decays slower than the bounded distribution, this probability for n = 3 
is found to be smaller than that for the uniform case. Thus the ‘reality’ is determined by both the tail and near-zero 
behavior of the probability distribution from which the matrix elements have been chosen. For the broad class of 
smooth distributions we have considered, we did not find the probability of real eigenvalues to lie outside the range 
[5/8, 7/8], which appear to be natural boundaries for sufficiently smooth distributions. 

Numerical results discussed above also shows that this remains valid even after taking a product of random matrices. 
It is quite surprising that the curves for probabilities of the different distributions do not cross each other. An intuitive 
explanation of this is yet to be figured out. We also measured the distribution of the matrix elements of the product 
matrix, and find that it decays slower (but faster than any power law) as the number K in the product increases. This 
observation is consistent with the result for single matrices (with finite moments), namely, that faster the distribution 
decays, smaller is the probability that all eigenvalues are real. However this is only part of the explanation as the 
probability p^^ N always stays below unity when correlations between matrix elements are ignored, and in order to 
approach unity, correlations are found to be crucial. It maybe stated that despite extensive and exact results obtained 
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so far there is not much understanding about why the eigenvalues tend to become real under the usual matrix product, 
and more work in this direction is clearly desired. 
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Appendix A: Detailed derivations for symmetric beta distribution 


When /i = 0 and v = 2k is a positive integer, the convolution can be written in the form of a finite series: 


q(z) = 0(2- |z|) 


2fc + 1 


2 2k 

E 

r=0 


2k 


H*l) 


2k 


_ r l-(lz|-l) 2fc+r+1 
2k + r + 1 


(Al) 


Integrations can be done using the hyperbolic variables as the combination xy appears exclusively. To outline the 
method used, hyperbolic coordinates (v, w ) where x = ve w and y = ve~ w are useful. Thus the variable w does not 
appear in the transformed expression. The range of v is [0,1], while for a given v, the range of w is [Inn, — Inti]. The 
integration over w is easily carried out first, and noting that the Jacobian of the transformation is 2u, leads finally to: 


2k 


P, I 2 V = 0, 1 / = 2fc) = 1 - (2fc + l) 4 x Y, 


r—0 

2/c+r+l 


2 k\ 2 2k ~ r+1 


(-l) r 


r J 2k + r + 1 \ (2k — r + 1)(6 k — r + 3) 2 


E 

;=o 


2k + r + 1 

l 


(- 2) 1 


(2k — r + l + l)(6/c — r + l + 3) 2 


(A2) 


As special cases: 


E 2 V = 0,u = 2)= = 0.631023, P^(p = 0,v = 4)= E” « 0.628361. 


72144072 


(A3) 


The case of odd integer v is also accessible, however we only state the result when v = 1. The convolution is now 
piecewise continuous on \z\ < 1 and 1 < \z\ < 2, and the probability of real eigenvalues is 


_(i), n ,, 3653 In 2 

2,2 (M- - v ~ ) - 576Q + 240 


0.63709. 


(A4) 


To understand the behavior of the probability of real eigenvalues for negative v, it is convenient to write the 
convolution for the distribution Po >v (x) defined by Eq. m as 


/ +00 

p(x)p(z — x) dx 

-OO 


v + 1 


r+l 


(A5) 


\ u \z — x\ v Q(l — \z — x\)dx. 


Note that \x\ < 1 and \z — x\ < 1 implies that \z\ < 2. From the symmetry q(z) = q(—z) it suffices to consider 
0 < z < 2. Using the hyperbolic coordinates, x = ve w and y = ve ~ w , we get 


n 1 /*— In v p 2 v 

Pp 2 = 1 — (v + l) 2 / 2 vdw 2v / dw dzq(z) 

J 0 J In v J 0 

/»1 p2v pi 

= l + (l + vf v 1+2v \nv dz \x\ v \x-z\ v Q(l-\x-z\)dxdv. 

Jo Jo J -1 


(A6) 


Except when v is an even integer, some care must be taken as q(z) is piecewise continuous in the two intervals [0,1] 
and [1,2]. For v = —1/2, we have 


pW (p = 0, v = -1/2) = 4(41 - 7T - 2 In 2) « 0.759836. 
’ 48 


(A7) 


To find the probability as v —> —1, consider the two inner integrals over z and x in Eq. (1A6I) . On considering 
various cases arising due to the Heaviside theta function and performing the integral over z first, we get 


_ n 3 r(i)T(2 + «/) l + 4(l + i/)ln2 


PM(p = 0 , 1 /) = -- 


F(| + v ) 


42(2+^) 


— J\ + J 2 , 


where 


T 1 = ( l + ,) 2 (J dv v^^lnv (l-2u) 1+1/ -^ dv v 1+2u lnu (2v - l) 1+,y ^ 


(AS) 


(A9) 
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and 


/ ,1/2 f l- 2 v ,1 ,i \ 

J 2 = (1 + v) 3 I 2 / dv u 2+2iy I 11 v / dx (a:(2u + x))" + / dv v 1+2l ‘ lnu / dx x v (2v — x) 1+ " ) . (AlO) 
y do do dl/2 J 2 v -1 J 


It is clear that when 1 / —> — 1, the second term on the right hand side of both J\ and J 2 is zero. The first integral in J\ 
contributes —1/4 when the integral is carried out by expanding the integrand around v = 0. However, the dominant 
contribution to the first integral in J 2 comes when both x and v are small. Then it is useful to split the integral over 
v from 0 to 1/4 and 1 /4 to 1/2 and it turns out that the former integral with x lying in the range 0 to 2v contributes 
— 1/16 to J 2 . Using these results and setting v = — 1 in the second term in Eq. (jA8j) . we finally obtain the desired 
result, viz., (/i = 0, v —> —1) = 7/8. 


Appendix B: Detailed derivations for smooth bounded distributions 


For the case rj = 1 (parabolic distribution), the convolution is still calculable as: 


q{z) = —(2 - M) 3 (4 + 6 |*| + z 2 )Q( 2 - |s|). 


(Bl) 


The integrals in Eq. ([5]) are then elementary and are easily done with mathematical packages, and yield the probability 


P$(V = 1) 


489341 

705600 


0.69351. 


(B2) 


The case for 77 = 2 leads to longer expressions, but with Mathematica, it is easy to evaluate the convolution and the 
resultant probability as 


9{z) = ^(2 - M) 5 (16 + 40|*| + 36z 2 + 10 |z | 3 + z 4 ) 0(2 - |*|), 


(B3) 


and 


pi% = 2 ) 


180521487191 

258564354048 


0.698168. 


(B4) 


Similar analysis yields the probability for 77 = 5,10,20,50 as 0.702769,0.704785,0.705906,0.706616 respectively. A 
monotonic convergence of the probability to l /\/2 as 77 increases is then a reasonable indication from these calculations. 


Appendix C: Detailed derivations for Gaussian distribution 


For the Gaussian distribution (1231) . the convolution with itself is also a Gaussian distribution but with variance of 
2 rather than 1: q(z) = e~ z /, 4 /( 2 y / 7 r). Application of Eq. ( 0 , along with the usage of hyperbolic coordinates results 
in 


p U) — 1 _ 
r 2,2 — 1 


poo p2^/xy 

[ e^ 2 ^ 2 )/ 2 f e- z2 /*d, 

Jo Jo 


r 3/2 


dx dy 


(Cl) 


1 poo POO p2v 

= 1 --oTo 2vdv e~ y2cosh2w dw / e-^^dz. (C2) 

7T J / 2 Jo J_ 00 _ Jo 

The integral over v becomes elementary when the burden of a finite range of integration is shifted from the * variable 
to the variable v 2 . That is, the v integral is done first following writing the * integral as e~ z / 4 0(—* + 2v) over the 
interval [0,oo). The probability of real eigenvalues simplifies to 


P w = i_j_ r dw _ 1 . 

2,2 nV2Jo cosh 2w cosh w y/2 


(C3) 


This follows as the integral can be evaluated as 7 t(\/ 2 — l)/2 either with packages such as Mathematica, or with the 
use of the residue theorem. This is applied to a rectangular contour with base on the interval [—L,L\, and the top 
on the interval with Im(ui) = iir and Re(w) € [—L,L\. This region encloses three poles of order 1 at 77 t/4, in/2 and 
i3n/4. In the limit of large A, the contour integral is twice what is required. 
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If we follow the path of the characteristic function, we have p(u>) = e u . Computing the integral in Eq. © yields: 

2 


pW = { — 

2 ’ 2 ^ o n\(2 n +iy 


2 (-1)^ 


r|2 + ? 

' 2 4 


(C4) 


A numerical approximation to Eq. (IC4I) using Mathematica yields a value which is exactly equal to l/\/2. Also, using 
the properties of Gamma functions, it is possible to modify this series to obtain 


ID i - ( -i r r(n + i)r(f + l) 

2.2 27T 2-^ n\ 

n—0 


r (f+ !) 


which has already been shown to be exactly equal to l /\/2 in El- 


Appendix D: Detailed derivations for Gamma distribution 


For 7 = 1, we have the so-called Laplace distribution p{x) = (l/2)e~\ x ' for which p 7 (0) is nonzero. Using the 
linear transformation formula (15.3.3) of [23| and the definition of the Gauss hypergeometric function 2 ^ 1 ( 0 , &; c; z) 
in Eq. ©, after some algebra, we obtain 

— D — 1 - - / duj(l + uj 2 ^- 7/2 
/ 0 


Pr 


1 r °° 

;,2 (7 =i) = 1 - ^ J o duj 0 + w2 )~ 


= — « 0.733333. 
15 


(Dl) 

(D2) 


For integer 7 , similar steps can be used to find an analytical expression for ■ 
The convolution is given by 


q(z) = 


4 r ( 7 ) 2 J-a 


■ - zp-HxP-'e-Me-^-^dx. 


(D3) 


Note that the general symmetry q(z) = q{—z) holds, and we can therefore restrict attention to z > 0. This simplifies 
the computation somewhat and leads to Thus the convolution is itself the sum of a Gamma distribution and a Bessel 
distribution [28]. If the distribution was not symmetrized, that is the distribution were the usual Gamma distribution 
on [0, 00 ), only a term proportional to the first will be in the convolution. The second term represents the convolution 
of the distributions on the opposite sides of the maxima. Now to evaluate the probability of real eigenvalues of a 2 x 2 
matrix with entries from the symmetric Gamma distribution, we use the first equality in Eq. ©, and obtain 


1 1 pOO /* OO pOO 

2,2 ( 7 ) = -Z + 7^2 / dx dyixyy^e^-y dzq(z). 

* -L w) Jo Jo Jz^/xy 


(D4) 


Once again, using the hyperbolic coordinates, x = ve w /2,y = ve w /2 and performing the integral over w results in 
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p 00 poo 

■ / dv u 2 "'^ 1 Kq{v) / dzq(z). 
Jo Jv 


2 227-2r(y)2 J Q 

We first consider the parameter regime 7 < 1. Apart from 7 = 1 the other “easy” case is 7 = 1/2 where we find 


(D5) 


^,2 (7 = 1/2) = ^ + l « 0.784155. 


(D 6 ) 


In fact, evaluation of such integrals with packages like Mathematica, even in this case, return unevaluated hypergeo¬ 
metric functions. However, in this case, one can simplify the integrals directly (hypergeometric functions encountered 
do not seem to have known identities, at least to our knowledge). For instance, in the evaluation with the Bessel 
distribution part of the convolution, one needs 
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/•OO /»00 -t /»C 

J J K 0 (v)K 0 (v + z) dv dz = J 


udu 1 


sinh u 8 


(D7) 


which is obtained by using the integral representation of the Bessel function of zeroth order and simplifying. The 
integral with the first term of q(z) is more easily evaluated to 1 /( 2 - 7 t). Thus the probability is 1/2 + 1/(27 t) + 1/8 as 
stated above. Numerical evaluation of these integrals is also tricky as 7 decreases toward zero. 
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For 7 > 1, the distribution = (1/2) + ^ 1 ( 7 ) + ^ 2 ( 7 ) where 


a m = ^ 1 r(fc + 2 7 ) 2 3 _ r(7 + ^) 2 1 r(fc + 7 ) 2 

2 4 5 6 tf( 7 ) 2 ^ 2^!r(fc + 2 7 +i)’ 2K1> 4v^r( 7 ) 2 ^ fc! r(fc + 27 + i) ■ 


(D 8 ) 


The sums above maybe written in terms of hypergeometric functions, but they do not appear to be simple expressions. 
To state some other exact values for the probability of real eigenvalues: 


pi 3(2) 


10259 

15015’ 


^$(3) 


640561 

969969' 


The sum A\ ( 7 ) for large 7 can be approximated by 


a x ( 7) w V -L_I_ (27 + fc)! 

UfJ 2 4 t'( 7 !) 2 ^2"v (2 7 + kf/ 2 k\ 


(D9) 
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where we have used that the ratio x\/{x + 4)! ~ x 172 for large x. Employing the Stirling’s approximation s! ss 
V / 27rs(s/e) s in the above equation, and replacing the sum by an integral, we obtain 


^ 1 ( 7 ) 


1 f 2 " 1 dk 27 

W^Jo + fc exp 



(Dll) 


where 


<7(2;) = ln(l + x) + a;ln(l + x x ) — (1 + x) In 2. 


(D 12 ) 


On evaluating the integral in Eq. (1D11I) by saddle point method, we obtain ^ 1 ( 7 ) = (l/ 8 )erf(-y/ 7 / 2 ) which, for large 
7 , gives 1/8. An analysis similar to above also shows the leading order correction in Ai(y) to be (16 v / 27T7) _1 , and 
41 . 2 ( 7 ) 1 ° be exponentially small in 7 . 


Appendix E: Distribution of the product for the probability distribution Eq. (1151) 

For the distribution Po. v {x) defined by Eq. m, the distribution of the product is given by 

/1 + i/\ K f 1 f 1 // 

Pk{zk) = J J dxi ■ J dxK \xi\ v -\x K \ v S(z K - X\xi) 


1 - f 1 f 1 K 

= t;( 1 + v ) K / dx 1 - / dx K x"...x v K 8{\z K \-T\xi), 
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where the last equation follows on noting that one half of the 2 K cases corresponding to the sign of the set { 27 } 
contribute equally to positive zk- On carrying out the integral over xk , we have 
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Pk(zk)=- (^ + ir \zk\ v / —W , — - 
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\ x k\ XK-1 
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~A 


, \zk\ < 1- 


(E3) 


(E4) 


The above result can be obtained by either transforming the problem of the distribution of product to that of sum 
by writing In 27 = yt, or directly using Mellin transforms j26j. 
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